Swarming and swirling in self-propelled polar granular rods 
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Using experiments with anisotropic vibrated rods and quasi-2D numerical simulations, we show 
that shape plays an important role in the collective dynamics of self-propelled (SP) particles. We 
demonstrate that SP rods exhibit local ordering, aggregation at the side walls, and clustering absent 
in round SP particles. Furthermore, we find that at sufficiently strong excitation SP rods engage in 
a persistent swirling motion in which the velocity is strongly correlated with particle orientation. 

PACS numbers: 45.70.Qj, 05.65.-l-b 



Large-scale structures emerge spontaneously in sys- 
tems of interacting SP biological objects such as flocks 
of birds, schools of fish, amoebae colonies, as well as in 
multi-robot swarms [J, Chemotaxis and field gra- 
dients can lead to non-equilibrium aggregation_|3| , and 
hydrodynamic interactions can cause vortices Such 
observations prompted a discrete-time, discrete-element 
model [sl where SP point particles { "boids" ) align their 
velocities with the average velocity of other particles 
within a certain fixed-size neighborhood. This model 
predicts a spontaneous phase transition to a long-range 
ordered state as the noise (temperature) of the system is 
reduced below a critical value, however the exact nature 
of the transition is still a matter of debate @ . Contin- 
uum hydrodynamic-type field models for a population 
of SP particles have been derived either general symme- 
try arguments Q or directly from microscopic interac- 
tion rules These models allowed for detailed pre- 
dictions of the correlation properties within the ordered 
state. However, both these models did not explicitly 
take into account the finite size and shape of interact- 
ing particles. On the other hand, there have been rapid 
advances in the theory of "active nematics", or popula- 
tions of inelastic ally interacting rods, both polar [§, [l3| 
and apolar [III, [13 ■ These models predict onset of a 
nematic order when the coupling strength of particle 
density becomes sufficiently high, furthermore, cluster- 
ing of apolar rods can lead to giant density fluctuations. 
Clustering of polar rods was recently found in numer- 
ical simulations On the experimental side, there 
has also been growing interest in the nonequilibrium dy- 
namics of driven granular rods. Symmetric rods in a 
vibrated container have been shown to form nematic or 
tetratic order and under certain conditions exhibit per- 
sistent swirling [13, [3, [i3| , and giant number fluctua- 
tions . At higher density, rods begin to bounce on one 
end and travel in the direction of their tilt due to friction 
at the contact between the rod and the substrate 17 1. 
Collectively, these rods spontaneously form large scale 
vortices El, [11. 



Here, we study experimentally and numerically the col- 
lective dynamics of vibrated polar granular rods which 



FIG. 1: The distributions of displacements in a time interval 
r (a) parallel, and (b) perpendicular to axis of the polar rod 
(F = 2). Inset: Schematic of the polar rod. Arrow shows the 
direction of net motion. 



have a non-symmetrical mass distribution but retain 
their shape symmetry. Such a rod on a vibrated surface 
moves towards the lighter end. When many such rods 
are placed inside a vibrated container, for weak excita- 
tion they aggregate over time at the boundaries. When 
the magnitude of excitation is increased, aggregation at 
the boundaries is reduced, and coherent structures are 
found in the bulk of the container. In particular, swirls 
can be identified in time averaged velocity fields, the flow 
and the rods are aligned, and signatures of incipient clus- 
tering can be observed. To augment these results and 
extend them toward larger system sizes, we perform nu- 
merical simulations using a discrete-element molecular 
dynamics algorithm. In particular, we show the impor- 
tance of particle aspect ratio and driving fluctuations on 
the observed pattern formation. 

About 10"^ polar rods were built using white hollow ny- 
lon cylinders of length / — 9.5mm and diameter d = 4.76 
mm, so the aspect ratio of the rods was fixed at 2 (see 
Fig. [2 a). Inset). Solid steel cylinders of length 4.75 mm 
and diameter 2.5 mm were placed snugly in one end of 
the nylon tube, which resulted in the center of mass be- 
ing displaced by O.ll from the geometrical center of the 
rod. The total mass of the assembly was 2.20 x 10^"* kg. 
The steel inserts also made the corresponding ends to 
appear somewhat darker and were used to identify the 
"polarity" of the rods. The particles were placed on a 
flat anodized aluminum container of radius R — 30d. 
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ing the arguments developed for symmetric rods and 
dimers 
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FIG. 2: (a) Rods migrate and aggregate at the boundaries of 
a container for modest excitations (A'^ = 500, F = 2). (b) Ag- 
gregation reduces and a homogeneous distribution is observed 
as excitation is increased. (A*' — 500, F = 4). (c) Area fraction 
p{r/R) as a function of distance r from the center of the con- 
tainer with radius R for A*' — 900 averaged over 100 frames at 
10 frames per second, open symbols: the results of numerical 
simulations for the same system parameters; (d) Simulations 
show decrease of clustering as Ar is reduced (F = 2) . 



The container was vertically vibrated using an electro- 
magnetic shaker with sinusoidal waveform at frequency 
75 Hz and varied driving acceleration T (scaled by the 
gravity acceleration) from to 5. A digital camera with 
the spatial resolution of 1000 x 1000 pixels was used to 
image the motion of rods inside the container. 

First, we studied the motion of a single rod bouncing 
on the vibrated plate away from side walls. For F > 1.5, 
the rod shows a robust net motion in the direction of the 
lighter end of the rod while taking some apparently ran- 
dom steps in the other directions as well. A movie of the 
typical motion is contained in the Supplementary mate- 
rial. By cross-correlating the intensity distribution of the 
image of the rod with a mask, we automated finding the 
position and the orientation (measured by the angle (f) to 
a fixed axis) of the rod in each frame. By measuring the 
change in position over time interval r, the magnitude 
of the rod velocity v, and its direction with respect 
to a fixed reference were obtained. The probability dis- 
tribution functions (PDF) for the displacement parallel 
to the rod vt cos{9 — (j)) and perpendicular to the rod 
VTsm{9 — (f>) are plotted in Fig. [Ha,b) with 3 x 10^ sets 
of measurements. While the PDF in the perpendicular 
direction are centered at zero, the broader PDF in the 
parallel direction are clearly shifted from zero, and this 
shift grows as r is increased. The mean and the rms ve- 
locity increase with F in our system (see Supplementary 
materials.) By imaging from the side, we find that rods 
undergo short collisions with the bottom of the container 
once every few cycles at random phases of the cycle (see 
Supplementary materials.) 

The physical mechanism for the observed directed mo- 
tion in our polar rods can be understood by extend- 



l20j . During a typical collision of a particle 
with a horizontal plate, a large but short impulse of fric- 
tional force at the contact point causes horizontal parti- 
cle displacement after the collision. When a symmetrical 
(apolar) rod bounces symmetrically on a vibrated plate, 
the net displacement after many collisions is absent, but 
for an asymmetric mode of vibration (as in Ref. |20j|) 
or for an asymmetric particle (as in the present study), 
there is a non-zero net horizontal motion. In the case 
of polar rods, since the center of mass is displaced from 
the geometrical center, the heavy end collides more often 
with the plate, and the rod on average travels in the di- 
rection of the light end. It can be shown that the average 
horizontal velocity of the rod translation is proportional 
to the amplitude of the vertical speed of the container, 
and indeed we observe that the mean velocity increases 
with F. 

The collective motion of polar rods was studied by 
placing the rods randomly initially inside the container 
and then vibrating with various F. (Movies included in 
Supplementary materials.) For low F ~ 2, rods were ob- 
served to migrate to the boundary of the container and 
aggregate in about 30 seconds. An example is shown in 
Fig. [2Ia). Not all rods aggregate at the boundaries, as 
some rods gradually rotate and escape from the dense 
cluster at the boundary back into the middle of the 
container. As F is increased, so do fluctuations, and 
the aggregation at the boundaries becomes less and less 
pronounced. Although spatio-temporal density inhomo- 
geneities persist, the time- averaged number density of the 
polar rods appear more or less uniform across the cavity 
for F > 3 (see Fig.lKb)). 

Next, we performed "molecular dynamics" simulations 
of polar rod motion and interaction. We did not simu- 
late the details of the vibrational transport of bouncing 
rods, but instead assume that the rods were confined to a 
horizontal plane. A force acts on each rod along its (hori- 
zontal) axis in the direction of the lighter end. This force 
was assumed to be random, with a mean F and variance 
V. In addition to the driving force, we assumed that rods 
experience velocity-dependent friction with the substrate 
and inelastic collisions with other rods. F and V were 
tuned so the displacement distribution for a single rod fits 
the experimental data for a given F. In the numerics, the 
rods had a form of spherocylinders, which helped in mod- 
eling contact forces. The interaction forces among rods 
were calculated via the interaction between viscoelas- 
tic virtual spheres of diameter d centered at the closest 
points between the axes of the spherocylinders [l3| . Nor- 
mal forces were computed using Hertzian spring-dashpot 
model, and dynamic Coulomb friction was assumed for 
tangential forces. We did not add random forcing in the 
direction perpendicular to the axis of the rod, so the rods 
could only change their direction by colliding with the 
walls or other rods. 
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FIG. 3: (a) The standard deviation of the number of rods An 
versus mean number of rods n inside a circular area at the 
center of the container (F = 3). Open symbols correspond to 
the simulations for a larger system size Rl ~ 2.5Rs but the 
same density. We consistently observe that An grows faster 
than y^n, which is a signature of a clustering regime. Local 
orientational order Q for nearest neighbors from simulation 
data as a function of (b) density, and (c) aspect ratio. 



We first performed simulations which matched the ex- 
periment both in terms of number of rods and the sys- 
tem geometry. We used aspect ratio Ar = 2 and im- 
posed elastic boundary conditions on a circle of radius 
Rs = M.2d 0. For F = 0.25, V = 0.16 which corre- 
spond to r = 2, we find that rods tend to aggregate at 
the boundary in agreement with experiment. As F and V 
are increased, the aggregation at the boundaries dimin- 
ishes, also in accord with the experimental observations 
(see numerical movies in the Supplementary material). 
To illustrate and compare the aggregation of rods in the 
experiments and simulations, we plot the projected rod 
area fraction p as a function of distance from the center 
r in Fig. (He). 

Clustering at the walls is not simply the consequence 
of inelastic collisions. Indeed, aggregation doesn't occur 
at the boundary for small Ar = 1.1 (Fig. El^d)), which 
indicates that particle shape affects aggregation. When 
fluctuations are small (at small T), rods have a much 
lower probability of turning around and leaving the wall 
than spherical particles, and so they are trapped near the 
wall for a long time. 

In order to characterize the density fluctuations inside 
the container, we obtain the standard deviation An and 
the mean n of the number of rods in areas of different 
sizes by averaging over many realizations (see Fig.[3l[a)). 
The distributions were obtained by averaging over 1500 
frames after the system reached a statistically stationary 
regime, and we restricted the area of measurements to 
r/R < 0.7 to minimize boundary effects. The data is 
systematically higher than i/n. In fact, they are better 
described by the slope 7/12 which is predicted by the dy- 
namic XY model [7| in the nematic state. At very high 



values of n the standard deviation drops down, as should 
be expected since the number of rods becomes compa- 
rable with the total number of rods in the container. 
In the numerical simulations which are also plotted in 
Fig. Efa), we examine a larger system with Rl = 2.5Rs 
and almost an order of magnitude greater number of 
rods. The deviations from ^y-n are also clearly present 
in this larger system. It is interesting to contrast these 
results with "giant" (An ^ n) ffuctuations reported for 
apolar rods [Til, [3] . Although rods in our system have 
apolar shape, they have mass anisotropy which renders 
them polar and self-propelled under external vibration. 
This polarity appears to destroy the emergence of giant 
density fluctuations in agreement with earlier theoretical 
work [3]. 

Although global orientational order is clearly absent 
in our system, there is a strong evidence of the local 
orientational order at sufficiently high density of rods. 
We can characterize this ordering by computing a local 
orientational order parameter Q which we define as Q = 
(cos 20) where & is the angle between directors of a rod 
and its nearest neighbor and brackets indicate averaging 
over all the rods in the container and time (see Fig.[3jDc). 
Parameter Q is similar to the local orientational order 
parameter S introduced in Q and shows significant local 
orientational order present in our system at high enough 
p and Ar- 

Collective motion of rods in the container is masked 
by the strong random fluctuations, especially at high F. 
To reduce these fluctuations, we divided the held of view 
into 2d X 2d boxes and averaged the velocity fleld over the 
box area and over a r = 5 second time interval. An ex- 
ample of the obtained velocity fleld is shown in Fig.jH^a). 
This procedure reveals a number of streams and swirls. 
Numerical simulations for similar parameters also show 
swirl-like structures [see Fig. |4ljb)]. The coherent struc- 
tures become more pronounced when the system size is 
increased [see Fig.Hl^c)]. These structures are reminiscent 
of swirls obtained with apolar particles driven by the sub- 
strate [13, 15 1 . However this is not entirely unexpected 
since the tensor order parameter for apolar particles in 
2D has only two independent components and the corre- 
sponding order parameter equation can be reduced to a 
pseudo- vector form jl5| which is similar to a vector order 
parameter equation for polar systems 13, 3 21 1 . 

To quantify the structure of swirls, we plot in Fig.|4jd) 
the spatial velocity correlation function C{r) = (vi • 
V2)/(|vi I |v2|) for a rod with velocity Vi and a rod with 
velocity V2 separated by distance r. The correlations de- 
cay over a distance of a few rod lengths which conflrms 
the lack of the long-range order in the system. However, 
the structure of the velocity field is strongly correlated 
with the orientation of the rods. We computed the distri- 
bution of the angle between the direction of the velocity 
field 9 in and the mean orientation within a {2d x 2d) 
box both in experiment and numerical simulations, see 
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rected motion in realistic self-propelled rods which affects 
the phenomenology of their collective dynamics. 
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FIG. 4: Swirling in time averaged velocity obtained by com- 
puting particle displacements after r = 5s: (a) experiment 
(r = 3iV = 900), (b) numerics (Fo = 1.0, p = 0.68), (c) 
Example of swirls observed in a larger numerical system 
(A^ = 5500, Rl ~ 2.5 Rs)- (d) Spatial velocity correlation 
function C(r) as a function of distance between two rods 
(r = 3, p = 0.31, p = 0.65); results of simulations are shown 
for small (.Rs = 34. 2d) and large Rl = 2.5Rs system sizes 
and for the same density p = 0.68. (e) The distribution of the 
angle between rod orientation and its velocity. 



Fig.mje). As seen in Fig.dl^e), there is a significant max- 
imum of this distribution at angle 0, which indicates that 
rods predominantly move along their axes. 

In summary, we have studied the collective dynam- 
ics of "self-propelled" polar rods with experiments and 
numerical simulations. The phenomenology differs qual- 
itatively from that of collective motion of both point-like 
self-propelled particles ^] (which show no tendency to 
aggregate near the walls and get involved in system-size 



collective motion) and apolar rods [ll|, Oj, |l6| (which 
exhibit giant density fluctuations). We observe aggrega- 
tion of rods at the boundaries because of the inability of 
rods to turn around and escape for high enough density 
under low noise conditions. As vibration strength and 
thus noise is increased, the aggregation reduces and a 
uniformly distributed state displaying local orientational 
order and swirls are observed. We observe greater than 
^/n density fluctuations which are in a qualitative agree- 
ment with model Q, but this agreement should not be 
over-emphasized since the model is directly applicable to 
a nematic regime. In our opinion, the observed devia- 
tion from the i/n regime for interacting polar rods is not 
accounted for by existing models and deserves further 
study. In conclusion, our findings elucidate an important 
and interesting interplay between the shape and the di- 
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